/*
 *     This file is part of Mobile GPS Logger.
 * 
 *     Mobile GPS Logger is free software: you can redistribute it and/or
 *     modify it under the terms of the GNU General Public License as
 *     published by the Free Software Foundation, either version 3 of the
 *     License, or (at your option) any later version.
 * 
 *     Mobile GPS Logger is distributed in the hope that it will be useful,
 *     but WITHOUT ANY WARRANTY; without even the implied warranty of
 *     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *     GNU General Public License v3 for more details.
 * 
 *     You should have received a copy of the GNU General Public License v3
 *     along with Mobile GPS Logger. If not, see <http://www.gnu.org/licenses/>.
*/

package cx.ath.skyflyer.gpslog.util;

/**
 * This class contains implementations of math functions that are missing from
 * the java.lang.Math in J2ME.
 * @author Ahti Legonkov
 */
public final class M {

    /**
     * This function does some magic to calculate the function and thanks to
     * that magic it's a bit more inaccurate but a lot faster than the other
     * implementation in pow.
     * @return a raised to the power of b.
     * @author don't know.
     */
    public static double pow_fast(final double a, final double b) {
        final int x = (int) (Double.doubleToLongBits(a) >> 32);
        final int y = (int) (b * (x - 1072632447) + 1072632447);
        return Double.longBitsToDouble(((long) y) << 32);
    }

    /**
     * @return a raised to the power of b.
     */
    public static double pow(final double a, final double b) {
        return exp(b*ln(a));
    }

    private static double ln_impl(final double x, final int n) {

        // 170 was determined with trial and error to be the limit
        // that yielded accuracy to 12 digits.
        if (n > 170)
            return 1.0;

        final double m = 1.0 + (double)((n-1)/2);
        return (double)n + m*m*x/ln_impl(x, n+1);
    }

    /**
     * @param y
     * @return natural logarithm of y.
     */
    public static double ln(final double y) {
        final double x = y - 1.0;
        return x / ln_impl(x, 1);
    }

    /**
     * @return e raised to the power of x.
     */
    public static double exp(double x) {

        //
        //   lim    pow(1+x/n, n) = pow(e, x)
        //  n->inf

        double n = (double)(1<<31);
        double res = 1.0 + x / n;

        for (int i = 0; i < 31; ++i) {

                res *= res;
        }

        return res;
    }

    public static double haversin(double x) {
        double sin = Math.sin(x/2);
        return sin*sin;
    }

    /** Arcus sine. */
    public static double asin(double x) {
        if (x < -1. || x > 1.) {
            return Double.NaN;
        }
        if (x == -1.) {
            return -PI_OVER_2;
        }
        if (x == 1) {
            return PI_OVER_2;
        }
        return atan(x / Math.sqrt(1 - x * x));
    }

    public static final double SQRT3 = 1.732050807568877294;
    public static final double PI_OVER_2 = Math.PI / 2;
    public static final double PI_OVER_6 = Math.PI / 6;
    public static final double PI_OVER_12 = Math.PI / 12;
    public static final double PI_OVER_180 = Math.PI / 180;

    /** Arcus tangent. */
    public static double atan(double x) {

        boolean signChange = false;
        // check up the sign change
        if (x < 0.) {
            x = -x;
            signChange = true;
        }

        // check up the invertation
        boolean Invert = false;
        if (x > 1.) {
            x = 1 / x;
            Invert = true;
        }

        // process shrinking the domain until x<PI/12
        double a;
        int sp = 0;
        while (x > PI_OVER_12) {
            sp++;
            a = x + SQRT3;
            a = 1 / a;
            x = x * SQRT3;
            x = x - 1;
            x = x * a;
        }

        // calculation core
        double x2 = x * x;
        a = x2 + 1.4087812;
        a = 0.55913709 / a;
        a = a + 0.60310579;
        a = a - (x2 * 0.05160454);
        a = a * x;

        // process until sp=0
        while (sp > 0) {
            a = a + PI_OVER_6;
            sp--;
        }

        // invertation took place
        if (Invert) {
            a = PI_OVER_2 - a;
        }

        // sign change took place
        if (signChange) {
            a = -a;
        }

        return a;
    }
}
